%
% Complete documentation on the extended LaTeX markup used for Insight
% documentation is available in ``Documenting Insight'', which is part
% of the standard documentation for Insight.  It may be found online
% at:
%
%     http://www.itk.org/

\documentclass{InsightArticle}

\usepackage[utf8]{inputenc}
\usepackage[dvips]{graphicx}
\usepackage{color}
\usepackage{minted}
\usepackage{float}
\definecolor{ltgray}{rgb}{0.93,0.93,0.93}
\usemintedstyle{emacs}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  hyperref should be the last package to be loaded.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[dvips,
bookmarks,
bookmarksopen,
backref,
colorlinks,linkcolor={blue},citecolor={blue},urlcolor={blue},
]{hyperref}


%  This is a template for Papers to the Insight Journal.
%  It is comparable to a technical report format.

% The title should be descriptive enough for people to be able to find
% the relevant document.
\title{Computing Textural Feature Maps for N-Dimensional images}

%
% NOTE: This is the last number of the "handle" URL that
% The Insight Journal assigns to your paper as part of the
% submission process. Please replace the number "1338" with
% the actual handle number that you get assigned.
%
\newcommand{\IJhandlerIDnumber}{3574}

% Increment the release number whenever significant changes are made.
% The author and/or editor can define 'significant' however they like.
\release{2.0.0}

% At minimum, give your name and an email address.  You can include a
% snail-mail address if you like.
\author{Jean-Baptiste Vimort$^{1}$, Matthew McCormick$^{1}$, François Budin$^{1}$ and Beatriz Paniagua$^{1}$}
\authoraddress{$^{1}$Kitware, Inc, Carrboro, NC}

\begin{document}

%
% Add hyperlink to the web location and license of the paper.
% The argument of this command is the handler identifier given
% by the Insight Journal to this paper.
%
\IJhandlefooter{\IJhandlerIDnumber}


\ifpdf
\else
   %
   % Commands for including Graphics when using latex
   %
   \DeclareGraphicsExtensions{.eps,.jpg,.gif,.tiff,.bmp,.png}
   \DeclareGraphicsRule{.jpg}{eps}{.jpg.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.gif}{eps}{.gif.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.tiff}{eps}{.tiff.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.bmp}{eps}{.bmp.bb}{`convert #1 eps:-}
   \DeclareGraphicsRule{.png}{eps}{.png.bb}{`convert #1 eps:-}
\fi


\maketitle


\ifhtml
\chapter*{Front Matter\label{front}}
\fi


% The abstract should be a paragraph or two long, and describe the
% scope of the document.
\begin{abstract}
\noindent
This document describes a new remote module implemented for the Insight Toolkit ITK
  \url{www.itk.org}, itkTextureFeatures. This module contains two texture analysis filters that are used to compute feature maps of N-Dimensional images using two well-known texture analysis methods. The two filters contained in this module are itkCoocurrenceTextureFeaturesImageFilter (which computes textural features based on intensity-based co-occurrence matrices in the image) and itkRunLengthTextureFeaturesImageFilter (which computes textural features based on equally valued intensity clusters of different sizes or run lengths in the image). The output of this module is a vector image of the same size than the input that contains a multidimensional vector in each pixel/voxel. Filters can be configured based in the locality of the textural features (neighborhood size), offset directions for co-ocurrence and run length computation, the number of bins for the intensity histograms, the intensity range or the range of run lengths. This paper is accompanied with the source code, input data, parameters and output data that we have used for validating the algorithm described in this paper. This adheres to the fundamental principle that scientific publications must facilitate reproducibility of the reported results.

\end{abstract}

\IJhandlenote{\IJhandlerIDnumber}
\newpage
\tableofcontents
\newpage
\section{Introduction}
\label{sec:intro}

Texture is an intuitive concept heuristically but difficult to define
precisely. It can be defined as series of homogeneous visual patterns that are
observed in certain kinds of materials. Humans describe texture through
qualitative concepts such as fine, coarse, granulated or smooth. These
descriptions are not precise and, in addition, they are not quantitative. The
texture quantification has been studied for a long time and many different texture quantification and analysis techniques exist. The work presented here explores texture quantification algorithms that provide a statistical description of the local texture of a 3D image and produces N-Dimensional texture maps. More importantly, the work here presented allows to obtain these outputs close to real-time.

ITK already contains tools for texture analysis such as
\doxygen{ScalarImageToTextureFeaturesFilter} or
\doxygen{ScalarImageToRunLengthFeaturesFilter}. Each of these filters can be
used in order to obtain a N-Dimensional texture characterization of an image.

However, those filters are implemented to describe images globally and they
are not optimized to be used to describe local texture features. Using the
existing filters to compute texture features iteratively for each pixels's
neighborhood, in order to build texture maps, quickly becomes too time
consuming. This is especially so when dealing with high resolution data, such
as micro computed tomography ($\mu$CT) or on large local neighborhoods of the input image.

The chosen solution, described in this article, consists in creating a new ITK
remote module (called itkTextureFeature) dedicated to the computation of feature maps
for N-Dimensional images. The filters implemented in itkTextureFeature
computes the exact same features as \doxygen{ScalarImageToTextureFeaturesFilter} and \doxygen{ScalarImageToRunLengthFeaturesFilter}. However, the new algorithms are optimized (particularly thanks to multithreading, \doxygen{NeighborhoodIterator}, \doxygen{ImageBoundaryFacesCalculator}) to be able to compute the feature maps much faster.

All the features available in itkTextureFeature are presented in
Section~\ref{sec:features}. Section~\ref{sec:filterUsage} describes the
filters specifications (templates, inputs, parameters) of each filter and how
to customize the use of these filters to each different texture analysis
application. Section~\ref{sec:examples} contain examples of code using
itkTextureFeature filters in Python and C++. Finally,
Sections~\ref{sec:results} and \ref{sec:conclusions} present several
scenarios, results and conclusions obtained with itkTextureFeatures.
\newpage
\section{Features Available}
\label{sec:features}

\subsection{Co-occurence features}
\label{sec:coocFeat}

The computation of co-occurrence features is based on the grey level co-occurrence matrix (GLCM) computed with \doxygen{itkCoocurrenceTextureFeaturesImageFilter} for each pixel's neighborhood. The GLCM matrix describes intensity organization of each pixels' neighborhood thanks to its second-order joint probability function \cite{coocFeat1,coocFeat2,coocFeat3,coocFeat4}. In the following computations the nomenclature is as follows:

\begin{math} g(i,j) \end{math} is the element in the cell \begin{math}(i,j)\end{math} of the normalized GLCM

\begin{math} \mu = \sum_{i,j}\nolimits i \cdot g(i, j) = \sum_{i,j}\nolimits j \cdot g(i, j) \end{math} is the weighted pixel average

\begin{math} \sigma = \sum_{i,j}\nolimits (i - \mu)^2 \cdot g(i, j) = \sum_{i,j}\nolimits (j- \mu)^2 \cdot g(i, j) \end{math} is the weighted pixel variance

\begin{math} \mu_t \end{math} and \begin{math} \sigma_t \end{math} are the mean and standard deviation of the row (or column, due to symmetry) sums

\doxygen{itkCoocurrenceTextureFeaturesImageFilter} computes the following features from the co-occurrence matrix:

\textbf{Energy} is a feature that measures the local uniformity of texture. The higher the energy value, the bigger the uniformity and organization of the texture.

\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_1 = \sum_{i,j}\nolimits g(i, j)^2
\end{equation}

\textbf{Entropy} is a feature that expresses the level of organization of a texture. A completely random distribution of grey-level intensities in the image volume would have very high entropy, while an image with the same grey-level across all pixels would have very low value of entropy.

\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_2 =
  \begin{cases}
     \sum_{i,j}\nolimits g(i, j)log_2g(i,j)  & \text{if } g(i, j) \neq 0 \\
     0                                       & \text{if } g(i, j) = 0
  \end{cases}
\end{equation}

\textbf{Correlation} is a feature that measures the linear dependency of gray level values in the co-occurrence matrix.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_3 = \sum_{i,j}\nolimits \frac{(i-\mu)(j-\mu)g(i, j)}{\sigma^2}
\end{equation}

\textbf{Inverse Difference Moment} (IDM) is a feature that measures the homogeneity of the image.  IDM will be low for in-homogeneous images, and a higher value for homogeneous images.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_4 = \sum_{i,j}\nolimits \frac{1}{1+(i-j)^2}g(i, j)
\end{equation}

\textbf{Inertia} or \textbf{contrast} is a feature that measures local grey-level variation in the GLCM matrix. If the neighboring pixels in the texture are very similar in their grey-level values then the contrast in the image is very low. Contrast is 0 for a constant image.

\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_5 = \sum_{i,j}\nolimits (i-j)^2g(i, j)
\end{equation}

\textbf{Cluster Shade} is a feature of the skewness of the matrix and is believed to be linked to perception of uniformity in the image. When this feature is high the image is asymmetric.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_6 = \sum_{i,j}\nolimits ((i - \mu)+(j-\mu))^3g(i, j)
\end{equation}

\textbf{Cluster Prominence} is a feature that is also related to the perceptual symmetry of the image. When the cluster prominence value is high, the image is less symmetric.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_7 = \sum_{i,j}\nolimits ((i - \mu)+(j-\mu))^4g(i, j)
\end{equation}

\textbf{Haralick's Correlation} is the original correlation measure designed by Haralick in 1973, and measures the linear dependence between pixels relative to each other.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_8 = \frac{\sum_{i,j}\nolimits (i,j)g(i,j)-\mu _t^2}{\sigma _t^2}
\end{equation}

\subsection{Run length features}
\label{sec:RLFeat}

The computation of the run length features is based on the grey level run length matrix (GLRLM) computed inside \doxygen{itkRunLengthTextureFeaturesImageFilter} for each pixel's neighborhood. A grey-level run is a set of consecutive, collinear picture points having the same grey-level value. The length of the run is the number of picture points in the run. The GLRLM matrix describes each neighborhood local texture. \cite{RLFeat1,RLFeat2,RLFeat3,RLFeat4}. In the following computations the nomenclature is as follows:

\begin{math} g(i,j) \end{math} is the element in the cell \begin{math}(i,j)\end{math} of the normalized GLRLM

\begin{math} i \end{math} is related the pixel intensity and \begin{math} j \end{math} to the length of the run

\doxygen{itkRunLengthTextureFeaturesImageFilter} computes the following features:

\textbf{Short run emphasis} (SRE) measures the distribution of short runs. SRE is expected large for fine textures.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_1 = \frac{\sum_{i,j}\nolimits \frac{g(i, j)}{j^2}}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Long run emphasis} (LRE) is a feature that measures distribution of long runs. LRE is expected large for coarse structural textures.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_2 = \frac{\sum_{i,j}\nolimits g(i, j)j^2}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Grey level non-uniformity} (GLN) measures the similarity of grey-level values through out the texture. The GLN is expected small if the grey-level values are alike throughout the whole texture.

\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_3 = \frac{\sum_{i}\nolimits (\sum_{j}\nolimits g(i, j))^2}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Run length non-uniformity} (RLN) is a feature that measures the similarity of the length of runs through out the image. The RLN is expected small if the run lengths are alike through out the image.

\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_4 = \frac{\sum_{j}\nolimits (\sum_{i}\nolimits g(i, j))^2}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Low grey level run emphasis} (LGRE) is orthogonal to SRE, and the value of the feature increases when the texture is dominated by many runs of low gray value.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_5 = \frac{\sum_{i,j}\nolimits \frac{g(i, j)}{i^2}}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{High grey level run emphasis} (HGRE) is orthogonal to LRE, and the metric increases when the texture is dominated by many runs of high gray value.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_6 = \frac{\sum_{i,j}\nolimits g(i, j)i^2}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Short run low grey level emphasis} (SRLGE) is a diagonal measurement that combines SRE and LGRE. The metric increases when the texture is dominated by many short runs of low gray value.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_7 = \frac{\sum_{i,j}\nolimits \frac{g(i, j)}{i^2j^2}}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Short run high grey level emphasis} (SRHGE) is orthogonal to SRLGE and LRHGE and increases when the texture is dominated by short runs with high intensity levels.

\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_8 = \frac{\sum_{i,j}\nolimits \frac{g(i, j)i^2}{j^2}}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Long run low grey level emphasis} (LRLGE) is complementary to SRHGE, it increases when the texture is dominated by long runs that have low gray levels.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_9 = \frac{\sum_{i,j}\nolimits \frac{g(i, j)j^2}{i^2}}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\textbf{Long run high grey level emphasis} (LRHGE) is the complementary metric to SRLGE and increases with a combination of long, high-gray value runs.
\begin{equation} \label{eqn:ShapeInfluenceTerm}
f_{10} = \frac{\sum_{i,j}\nolimits g(i, j)i^2j^2}{\sum_{i,j}\nolimits g(i, j)}
\end{equation}

\section{Filters usage}
\label{sec:filterUsage}

\subsection{itkCoocurrenceTextureFeaturesImageFilter}
\label{sec:coocFilter}

For each pixel of the input image, the itkCoocurrenceTextureFeaturesImageFilter will compute a serie of 8 co-occurrence texture features which will be contained in a vector. That way the output of the filter is a N-D image where each pixel will contain a vector of 8 scalars. Each texture map can be extracted from the output image afterward thanks to \doxygen{NthElementImageAdaptor}. By default the texture features are computed for each spatial direction and averaged afterward to provide rotationally invariant texture descriptors.

Template Parameters (if used in C++):
\begin{itemize}
 \item The input image type: a N-Dimensional image where the pixel type MUST be integer.
 \item The output image type: a N-Dimensional image where the pixel type MUST be a vector of floating points or an ImageVector.
\end{itemize}

Inputs and parameters:
\begin{itemize}
 \item An input image
 \item A mask defining the region over which texture features will be calculated. (Optional)
 \item The pixel value that defines the "inside" of the mask. (Optional, defaults to 1 if a mask is set.)
 \item The number of intensity bins. (Optional, defaults to 256.)
 \item The set of directions (offsets) to average across. (Optional, defaults to {(-1, 0), (-1, -1), (0, -1), (1, -1)} for 2D images and scales analogously for ND images.)
 \item The pixel intensity range over which the features will be calculated. (Optional, defaults to the full dynamic range of the pixel type.)
 \item The size of the neighborhood radius. (Optional, defaults to 2.)
\end{itemize}

\subsection{itkRunLengthTextureFeaturesImageFilter}
\label{sec:RLFilter}

For each pixel of the input image, the itkRunLengthTextureFeaturesImageFilter will compute a serie of 10 run length texture features which will be contain in a vector. That way the output of the filter is a N-D image where each pixel will contain a vector of 10 scalars. Each texture map can be extracted from the output image afterward thanks to \doxygen{NthElementImageAdaptor}. By default the texture features are computed for each spatial direction and averaged afterward.

Template Parameters:
\begin{itemize}
 \item The input image type: a N dimensional image where the pixel type MUST be integer.
 \item The output image type: a N dimensional image where the pixel type MUST be a vector of floating points or an ImageVector.
\end{itemize}

Inputs and parameters:
\begin{itemize}
 \item An input image
 \item A mask defining the region over which texture features will be calculated. (Optional)
 \item The pixel value that defines the "inside" of the mask. (Optional, defaults to 1 if a mask is set.)
 \item The number of intensity bins. (Optional, defaults to 256.)
 \item The set of directions (offsets) to average across. (Optional, defaults to {(-1, 0), (-1, -1), (0, -1), (1, -1)} for 2D images and scales analogously for ND images.)
 \item The pixel intensity range over which the features will be calculated. (Optional, defaults to the full dynamic range of the pixel type.)
 \item The distance range over which the features will be calculated. (Optional, defaults to the full dynamic range of double type.)
 \item The size of the neighborhood radius. (Optional, defaults to 2.)
\end{itemize}

\subsection{Recommendations}
\label{sec:recommendations}

Using the itkTextureFeature's filters with the default settings will lead, in all likelihood, to meaningless results. In addition, those results might be really time consuming to compute.

To obtain significant results, most of the parameters need to be carefully chosen depending on the scale of the texture and resolution of the input data, the other parameters and the significant information that need to be revealed by the output. For example the pixel intensity range should be adapted to the actual range of the input data in order to obtain a larger amplitude in the output feature maps. The radius (of the neighborhood) should be chosen depending on the size of the anomaly/object that needs to be detected in the input image. The distance range should correspond, in most cases, to the longer run length possible in the neighborhood (it depend of the spacing in the input image, the radius of the neighborhood and the offset directions).The choices of the number of bins should be adapted to the intensity/run length variation that want to be observed as well as the intensity/distance ranges specified. Finally the offset direction can be adapted if the anomaly/object that needs to be detected is specific to one known direction of the input image.

The usage of a Region Of Interest (ROI) mask is strongly advised, it will reduce the computation time by avoiding computing features for the noisy/background parts of the image.

In addition to the settings, particular attention should be payed to the input data. Please consider cropping the input are to contain only areas that will be interesting for the analysis. This will both help improving the computation time, thanks to a better distribution of the threaded regions and avoiding memory problems due to too large output data (considering that the output data is 8 or 10 times bigger than the input data).

The memory problem due to too large output data can also be solved by separating the output image containing all the feature maps into several images containing one feature map each thanks to the itk class \doxygen{VectorIndexSelectionCastImageFilter}.

\subsection{Python wheels}
\label{sec:PythonWheels}

Python wheels allow to easily install itkTextureFeatures and all its dependencies in order to have this texture filters ready to use in python code.They have been generated for the three main operating systems (Mac, Linux and Windows) and three versions of python (2.7, 3.5 and 3.6).

To install the python wheels use the following command-line: \$ pip install itk\_textureFextures

\newpage
\section{Practical examples}
\label{sec:examples}

\subsection{C++}
\label{sec:C++Ex}

\begin{minted}[baselinestretch=1,fontsize=\footnotesize,linenos=false]{cpp}
#include "itkRunLengthTextureFeaturesImageFilter.h"

#include "itkImage.h"
#include "itkVector.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkNeighborhood.h"

int main(int argc, char * argv[])
{
    if( argc != 10 )
    {
    std::cerr << "Usage: "<< std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputFilePath> <MaskFilePath> <OutputFilePath> ";
    std::cerr << " <NumberOfBinsPerAxis> <PixelValueMin> ";
    std::cerr << " <PixelValueMax>  <DistanceValueMin> ";
    std::cerr << " <DistanceValueMax> <NeighborhoodRadius> ";
    std::cerr << std::endl;
    return EXIT_FAILURE;
    }

    // Setup types
    using InputImageType = itk::Image< int, 3 >;
    using OutputImageType = itk::Image< itk::Vector< float, 10 > , 3 >;
    using readerType = itk::ImageFileReader< InputImageType >;
    using NeighborhoodType = itk::Neighborhood<typename InputImageType::PixelType,
      InputImageType::ImageDimension>;
    NeighborhoodType neighborhood;

    // Create and setup a reader
    readerType::Pointer reader = readerType::New();
    reader->SetFileName( argv[1] );

    // Create and setup a maskReader
    readerType::Pointer maskReader = readerType::New();
    maskReader->SetFileName( argv[2] );

    // Apply the filter
    using FilterType = itk::Statistics::RunLengthTextureFeaturesImageFilter
                                < InputImageType, OutputImageType >;
    FilterType::Pointer filter = FilterType::New();
    filter->SetInput(reader->GetOutput());
    filter->SetMaskImage(maskReader->GetOutput());
    filter->SetNumberOfBinsPerAxis(std::atoi(argv[4]));
    filter->SetPixelValueMinMax(std::atof(argv[5]),std::atof(argv[6]));
    filter->SetDistanceValueMinMax(std::atof(argv[7]),std::atof(argv[8]));
    neighborhood.SetRadius( std::atoi(argv[9]) );
    filter->SetNeighborhoodRadius(neighborhood.GetRadius());

    // Create and setup a writter
    using WriterType = itk::ImageFileWriter< OutputImageType  >;
    WriterType::Pointer writer = WriterType::New();
    std::string outputFilename = argv[3];
    writer->SetFileName(outputFilename.c_str());
    writer->SetInput(filter->GetOutput());
    writer->Update();

  return EXIT_SUCCESS;
}
\end{minted}
\normalsize

\subsection{Python}
\label{sec:pythonEx}

\begin{minted}[baselinestretch=1,fontsize=\footnotesize,linenos=false]{python}
import itk, sys

if len(sys.argv) != 7:
    print("Usage: " + sys.argv[0] + " <inputImagePath> <maskImagePath>"
                                    " <NumberOfBinsPerAxis> <PixelValueMin> "
                                    "<PixelValueMax> <NeighborhoodRadius>")
    sys.exit(1)

im = itk.imread(sys.argv[1])
mask = itk.imread(sys.argv[2])

filtr = itk.CoocurrenceTextureFeaturesImageFilter.New(im)
filtr.SetMaskImage(mask)
filtr.SetNumberOfBinsPerAxis(int(sys.argv[3]))
filtr.SetPixelValueMinMax(int(sys.argv[4]), int(sys.argv[5]))
filtr.SetNeighborhoodRadius([int(sys.argv[6]),int(sys.argv[6]),int(sys.argv[6])])

result = filtr.GetOutput()

itk.imwrite(result, "result.nrrd")
\end{minted}
\normalsize
\newpage
\section{Results}
\label{sec:results}

We present concrete use case scenarios of the different filters of itkTextureFilters in this section. We used itkTextureFilters to characterize subchondral bone texture in temporomandibular joint (TMJ) Osteoarthritis (OA). To date, there is no single sign, symptom, or test that can clearly diagnose early stages of TMJ OA. However, it has been observed that changes in the subchondral bone occur in very early stages of this disease involving subchondral bone structural changes (texture) in the subchondral bone (i.e. bone marrow).

The different tools presented in this document aid detecting and highlighting those texture variations in order to help clinicians to detect  TMJ OA earlier.

In the test case (figure \ref{fig:Scan}), the lower part of the condyle is healthy (normal bone trabeculae density) whereas the upper part is characteristic of a TMJ OA case (low bone trabeculae density).

\begin{figure}[H]
  \begin{center}
    \includegraphics[width=0.8\textwidth]{figures/Scan.eps}
    \itkcaption{CBCT of the test condyle: this condyle suffers of a lack of trabecula in the upper part}
    \label{fig:Scan}
  \end{center}
\end{figure}

\subsection{itkCoocurrenceTextureFeaturesImageFilter}
\label{coocResults}

The results exposed in this part were obtained by specifying the following parameters (the default parameters were used for the other ones):

\begin{itemize}
 \item Input data: Scan\textunderscore CBCT\textunderscore 13R.nrrd
 \item Input mask: SegmC\textunderscore CBCT\textunderscore 13R.nrrd
 \item Number of intensity bins: 10
 \item Pixel intensity range: [0; 4200]
 \item Neighborhood Radius: 6
\end{itemize}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/Energy.eps}
    \includegraphics[scale=0.3]{figures/Entropy.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's energy (left) and texture's entropy (right)}
    \label{fig:Energy&Entropy}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.25]{figures/Correlation.eps}
    \includegraphics[scale=0.25]{figures/InverseDifferenceMoment.eps}
    \includegraphics[scale=0.25]{figures/Inertia.eps}
    \includegraphics[scale=0.265]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's correlation (left), texture's inverse difference moment (center) and texture's inertia (right)}
    \label{fig:Correlation&Inertia}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/ClusterShade.eps}
    \includegraphics[scale=0.3]{figures/ClusterProminence.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's cluster shade (left) and texture's cluster prominence (right)}
    \label{fig:Clusters}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/HaralickCorrelation.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's Haralick correlation}
    \label{fig:Haralick}
  \end{center}
\end{figure}

It seems that, the energy (figure \ref{fig:Energy&Entropy}), the entropy (figure \ref{fig:Energy&Entropy}), the correlation (figure \ref{fig:Correlation&Inertia}), the inverse difference moment (figure \ref{fig:Correlation&Inertia}), and the inertia (figure \ref{fig:Correlation&Inertia}) are discriminating the healthy part of the bone trabecula from the unhealthy part of the bone trabecula. The detection of a significant variation in several of those feature maps could be interpreted in an early TMJ OA. In this particular case, the cluster shade feature map (figure \ref{fig:Clusters}), cluster prominence feature map (figure \ref{fig:Clusters}) and hahlick correlation feature map (figure \ref{fig:Haralick}) don't seem to discriminate different texture types within the TMJ condyle.



\subsection{itkRunLengthTextureFeaturesImageFilter}
\label{RLResults}

The results exposed in this part were obtained by specifying the following parameters (the default parameters were used for the other ones):

\begin{itemize}
 \item Input data: Scan\textunderscore CBCT\textunderscore 13R.nrrd
 \item Input mask: SegmC\textunderscore CBCT\textunderscore 13R.nrrd
 \item Number of intensity bins: 10
 \item Pixel intensity range: [0; 4200]
 \item Distance range: [0; 1.25]
 \item Neighborhood Radius: 4
\end{itemize}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/GreyLevelNonuniformity.eps}
    \includegraphics[scale=0.3]{figures/RunLengthNonuniformity.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's grey level non uniformity (left) and texture's run length non uniformity (right)}
    \label{fig:NonUniformity}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/ShortRunEmphasis.eps}
    \includegraphics[scale=0.3]{figures/LongRunEmphasis.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's short run emphasis (left) and texture's long run emphasis (right)}
    \label{fig:RunSize}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/LowGreyLevelRunEmphasis.eps}
    \includegraphics[scale=0.3]{figures/HighGreyLevelRunEmphasis.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's low grey level run emphasis (left) and texture's high grey level run emphasis (right)}
    \label{fig:GreyLevel}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/ShortRunLowGreyLevelEmphasis.eps}
    \includegraphics[scale=0.3]{figures/ShortRunHighGreyLevelEmphasis.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's short run low grey level emphasis (left) and texture's short run high grey level run emphasis (right)}
    \label{fig:ShortRunGreyLevel}
  \end{center}
\end{figure}

\begin{figure}[H]
  \begin{center}
    \includegraphics[scale=0.3]{figures/LongRunLowGreyLevelEmphasis.eps}
    \includegraphics[scale=0.3]{figures/LongRunHighGreyLevelEmphasis.eps}
    \includegraphics[scale=0.32]{figures/discreteFullRainbow.eps}
    \itkcaption{Texture's long run low grey level emphasis (left) and texture's long run high grey level emphasis (right)}
    \label{fig:LongRunGreyLevel}
  \end{center}
\end{figure}

Similarly to itkCoocurrenceTextureFeaturesImageFilter, some output feature maps of itkRunLengthTextureFeaturesImageFilter characterize bone texture and provide discriminant values of low and high density bone texture. Run length non uniformity (figure \ref{fig:NonUniformity}), short run emphasis (figure \ref{fig:RunSize}), long run emphasis (figure \ref{fig:RunSize}), high grey level run emphasis (figure \ref{fig:GreyLevel}), short run high grey level run emphasis (figure \ref{fig:ShortRunGreyLevel}) and long run low grey level emphasis (figure \ref{fig:LongRunGreyLevel}) seem to have discriminant properties and differentiate healthy and pathological bone trabecula. Here again some feature maps are not helpful in this particular case.

\newpage
\section{Conclusion}
\label{sec:conclusions}

This document presented a new fast and efficient tool to compute textural features in a  N-Dimensional image. These features are used to describe and detect variations in the image's texture. A lot of the described features are correlated to each other, so probably not using all of them at the same time is important. The variety of features available allow to detect a large type of variations, and using any combination of them might help texture characterization and discrimination. This method is currently used in a study aiming to create a new method to detect TMJ OA at early stages using subchondral bone texture as a biomarker.

\section*{Acknowledgements}

This work was supported by the National Institute of Health (NIH) National Institute for Dental and Craniofacial Research (NIDCR) grant R01EB021391 (Textural Biomarkers of Arthritis for the Subchondral Bone in the Temporomandibular Joint), NIDCR grant R01DE024450  (Quantification of 3D bony Changes in Temporomandibular Joint Osteoarthritis) and National Institute of Biomedical Imaging and Bioengineering (NIBIB) grant R01EB021391 (Shape Analysis Toolbox for Medical Image Computing Projects). 

We would like to thank Dr. Larry Wolford from the Baylor University Medical Center for kindly providing the bone specimens from which we obtained the scans used in the paper. We would like to thank Drs. Lucia Cevidanes, Erika Benavides and Antonio Ruellas at the University of Michigan School of Dentistry as well, for generating the CBCT scans that were processed with the filters presented in the paper.

We are also grateful for the support received from the ITK community.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Insert the bibliography using BibTeX
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\bibliographystyle{plain}
\bibliography{InsightJournal}


\end{document}

